BEYOND THE GROSS-PITAEVSKII EQUATION: 
GROUND STATE AND LOW ENERGY EXCITATIONS 

FOR TRAPPED BOSONS 

A. Polls 

Departament d'Estructura i Constituents de la Materia, 
Universitat de Barcelona 
Diagonal 647, Barcelona 08028, Spain 

A. Fabrocini 

Dipartimento di Fisica, Universita di Pisa 
and INFN, Sezione di Pisa, 1-56100 Pisa, Italy 



1. INTRODUCTION 

The recent experimental realization of Bose-Einstein condensation (BEC) of 
magnetically trapped alkali atoms has generated a huge amount of experimental 
and theoretical activity, with over twenty experimental groups in the world that can 
produce atomic condensates and more than one thousand articles on the subject 
pubhshed by now. The present status of the field has been recently reviewed by F. 
Dalfovo et al.[l]. 

The history of BEC gets back to 1924, when S. Bose derived the law for black- 
body radiation by considering the photons as a gas of identical particles. Immediately 
after, A. Einstein generalized Bose's ideas to an ideal Bose gas and predicted that the 
bosons would condensate in the lowest quantum state of the system at sufficiently low 
temperature. It took seventy years of experimental efforts and thechnological progress 
until the first atomic Bose condensate was achieved in 1995. In this year, groups at 
the University of Colorado and at the Massachusetts Institute of Technology (MIT) 
using a laser cooled and magnetically trapped dilute gas of ^^Rb [2] and ^^Na [3] 
atoms respectively, were able to demonstrate unambigously the occurrence of BEC. 

In those early experiments, the number of trapped atoms was relatively small 
(iV ~ 10'^ atoms) and the transition temperature was around 100 nK. The magnetic 
trap is well described by a harmonic oscillator potential, usually with cylindrical 
symmetry. However, througout this paper we will consider a spherical potential well 



confining the atoms. 

In order to have quantum effects, i.e., wave behaviour, we need a de Broghe 
wave length A = {2nh'^ /mT)^/'^ of the order of the distance between the atoms 
{pX^ ~ 1). On the other hand, the system should be kept dilute, therefore the 
critical temperature will be extremely low, of the order of nanokelvins. 

Up to now, the experimental conditions were such that the atomic gas was very 
dilute, i.e., the average distance between the atoms is much larger than the range 
of the interaction. As a consequence, the physics should be dominated by two-body 
collisions, generally well described in terms of the s-wave scattering length a. The 
case of a positive scattering length is equivalent to considering a very dilute system 
of hard spheres, whose diameter coincides with the scattering length itself. 

Typical scattering lengths are 53 A for ^^Rb and 28 A for ^^Na. On the 
other hand, the size of the trap is defined by the harmonic oscillator length oho = 
(h/mu)^^'^ which is of the order of 10^ A. The corresponding distance between the 
energy levels associated with this potential well is around 4 nK. For those initial ex- 
periments, a common ^^Rb atom density in the trap was p ~ 10-^^ — 10^^ atoms/cm^ 
giving an average inter-atom distance d ~ p~^l^ ~ lO'^A. Therefore, the effective 
atom size, defined by the scattering length is usually small compared to both the 
trap size and the inter-atom distance. The crucial parameter that defines the con- 
dition of diluteness is a; = pa^, which until very recently was kept rather small (i.e., 
X ~ 10~^). Under these conditions, the Gross-Pitaevskii equation [4], which assumes 
all the particles in the condensate, seems the logical tool to study those systems. 

The situation is somehow different in homogeneous liquid ^He. In this case, BEG 
manifests itself as a macroscopic occupation of the zero momentum state, measured 
by the condensate fraction, i.e., the fraction of the total number of particles in this 
state. However, there is only indirect evidence for this macroscopic occupation. 
Theoretical calculations and the analysis of inelastic neutron scattering data predict 
a condensate fraction of ~ 10 % [5]. Such a large depletion is an indication that '^He 
liquid is a strongly correlated system. 

There are two ways to bring x outside the regime of validity of the mean field 
description. The first consists in increasing the density and the second in changing 
the effective size of the atoms. Recent experiments have explored both possibilities. 
On one side they have reached a very high number of atoms in the condensate, ~ 10^, 
and on the other they have been able to change the scattering length of the atoms. 
This is the case of a recent experiment employing ^^Rb, where, by taking advantage 
of the presence of a Feshbach resonance at a magnetic field B ~ 155 Gauss, it was 
possible to vary the scattering length from negative to very high positive values. 
Under these conditions, effects beyond the mean field approximation are expected to 
be observable [6,7]. 

We will start by discussing a homogeneous system of Bose hard spheres and 
using the results to determine the regime of validity of the Gross-Pitaevskii equation. 
We will then present an extension of the GP equation [8], still in the framework of 
mean field theories, that allows for giving a first estimate of the expected corrections 
to the GP results in these new scenarios. Finally, we will briefly analyze the effects 
on the low collective excitations of the system. 



2. Homogeneous hard-sphere Bose gas 



We consider a system of N spinless bosons having mass m and described by the 
many-body Hamiltonian: 

I Kj 

The uniform system is studied in the thermodynamic hmit, N ^ oo and — > 00 
keeping the density, p = N/Q,, constant. The hard-spheres potential is defined as 
V{r < a) = 00 and V{r > a) = 0. 

Correlated Basis Functions (CBF) theory provides a very efficient way to handle 

the correlations induced by the interactions between the particles (for a review, see 
Ref. [9]). In its simplest version, one takes a Jastrow correlated wave function [10] 

*j(l,...,iV) = J]/(r.,), (2) 

i<j 

where the Jastrow correlation function, /(r), depends only on the interparticle dis- 
tance. Once the trial function is defined, the variational principle ensures that, if we 
are capable to calculate the expectation value of the Hamiltonian, 
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then i?cBF is an upper bound to the true ground state energy. The correlation 
function, /(r), is variationally determined by minimizing i?cBF- 

Eqbf may be calculated once the two-body distribution function, g{r), is known. 
In fact, the energy per particle can be written as 

e=^p J d^r g{r) 

In the particular case of hard spheres, the distribution function is strictly zero for 
r < a and the previous expression reduces just to the kinetic energy part 

^=-\pj d'r gir)^W^\nf{r). (5) 

g{r) may be evaluated by using cluster expansion and Hypernetted Chain (HNC) 
theory. HNC is an integral equation method which allows for massive summations 
of the cluster diagrams associated with g{r). 

The optimal choice for the Jastrow factor would be the one satisfying the Euler 
equation SEcBp/^f = 0. However, a less demanding and often effective approach 
consists in chosing a parametrized functional form of /(r) and in minimizing the 
energy with respect to the parameters. We adopt here the correlation function min- 
imizing the two-body cluster energy of a homogeneous Bose gas, with the healing 
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Figure 1. Energy per particle (in units of /i^/2ma^) for homogeneous hard spheres in 
function of x. The symbols correspond to the low-density expansion results obtained 
by keeping only the first term (LDq) or by adding the second (LDi) and the third 
(LD2) ones, and to the diffusion Monte Carlo (DMC) and HNC energies. The LD2 
energy at a; = 10~^ is E/N — 0.07 and lies outside the frame. 



conditions at a distance d {f{r > d) = 1 and f'{d) =0). d is taken as a variational 
parameter. For the hard-spheres case, /(r < a) — and /(r > a) = u{r)/r, where 
u{r) is the solution of the Schrodinger-like equation: —u" = K'^u. f{r) has the form 

[11] 



j.^^^ _ dsin[K{r - a)] 



(6) 



r sm[K(d — a)] ' 

and the healing conditions are satisfied through the relation: cot[K{d — a)] = {Kd)~^ 

An alternative calculation, based on perturbation theory in the expansion par 
rameter x — pa^, leads to the low-density expansion for the energy density[12]: 
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Up to these orders of the expansion, the details of the potential do not show up, 
and any potential with the same scattering length would give identical results. This 
universal behavior has recently been checked by a diflFusion Monte Carlo calcula- 
tion (DMC) [13], which provided the exact solution of the many-body Schrodinger 
equation. 

Fig. 1. shows the energy per particle in units of f? j2ma^ for homogeneous hard 
spheres as a function of x. The energies have been multiplied by lO^*-^'^-* at x = 



I0~^i~^'~^\ respectively. The figure compares the energies computed by retaining 
different terms in expansion (7). The LDq values contain the first term only, whereas 
LDi and LD2 are obtained by adding the second and third terms, respectively. The 
HNC results have been obtained disregarding the elementary diagrams (HNC/0) and 
using the correlation function of Eq.(6). The DMC results correspond to diffusion 
Monte Carlo calculations [13]. 

The agreement between the HNC/0 and the DMC results is excellent in most of 
the wide range of densities considered. However, there is a 5% disagreement at the 
highest x-value {x — 10~^). From this value on, the contribution of the elementary 
diagrams and a better optimization of the correlation should be probably taken into 
account. The LDq results are only accurate at very low densities, while LDi gives 
also a good representation of the exact DMC results. On the contrary, the addition 
of the logarithmic term spoils the agreement already at intermediate densities. 

In the next section, we will describe the trapped bosons by a local density ap- 
proximation (LDA) employing the homogeneous gas results. The local value of the 
parameter x in the trap will give an idea of the differences that we can expect by 
using the different energies reported in Fig. 1 as inputs to build the energy functional. 



2. Ground state of trapped hard spheres 



The energy functional associated with the Gross-Pitaevskii theory is simply ob- 
tained in the local-density approximation by keeping only the first term in the low 
density expansion (7): 
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where the wave function '^{r), in which all the atoms belong to the condensate, is 
normalized to N. By a functional variation of £?gp[*] one finds the Gross-Pitaevskii 
equation, 

/i^ _o moo 47r/i^a,^, , 
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*(r) = At*(r), 



(9) 



where ^ is the chemical potential. The CP equation has the form of a nonlinear 
stationary Schrodinger equation, and it has been solved for several types of traps 
using different numerical methods. 

The next logical step, in the spirit of LDA, is to include into the energy functional 
the next terms of the correlation energy expansion for the uniform system. According 
to the behavior of the different terms, shown in Fig. 1, it seems clear that it is 
reasonable to consider only the first correction, LDi. However, before proceeding 
further it is convenient to simplify the notation by expressing lengths and energies 
in harmonic oscillator (HO) units. The spatial coordinates, the energy, and the wave 
function are rescaled as r = auof, E = hujE, and *(r) = (iV/a^o)^^^*i(^)' where 
^'i(f) is normalized to unity. 

Using these new variables and taking into account the second term of the ex- 
pansion, we obtain the modified Gross-Pitaevskii (MGP) energy functional for the 



energy per atom, cmgp = Emgp/N, 
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and the corresponding modified Gross-Pitaevskii equation, 
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*i(r) =/^i*i(r), (11) 



where a = a/ano and \i is the chemical potential in HO units. 

In alternative, we have also used as local correlation energy the one provided 
by the uniform system HNC results. This option has the advantage that one is not 
limited to hard spheres, but, in principle, any type of potential for the two-body 
interaction can be considered. In this case, the local correlation energy V^^^, is given 

by 



yto?. = cJrpi(f)ete(pi). 



(12) 



where e^f^{pi) is the HNC homogeneous gas energy per particle at density pi. The 
of the energy gives the HNC correlated Hartree equation (CHhnc), 



dxioc 



^'i(r) = /ii^'i(r), (13) 



where we have introduced the scaled unities and the local gas parameter, x\oc{r) = 
pi(f)a3 = iVa3|*i(f)|2. 

The CP, MOP and CHhnc equations have been solved by the steepest descent 
method for an isotropic harmonic oscillator trap. 

Several relationships between the different contributions to the energy per atom 
or to the chemical potential exist, and are useful to check the accuracy of the numer- 
ical procedure. By direct integration of the CP equation, one finds 



A* = ekin + eno + 2e> 
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where Ckin = -|/ •i^r^'i V^^'i, eno = ^ / d^r\'^i\'^r^ and e\ll = f d^r27TaN\'^i\^ 
are the different terms contributing to the total energy. Further relationships can be 
obtained by means of the virial theorem. 



2ekin - 2eHO + 3e£J = 



(15) 



It is also important to notice that the dimensionless parameter characterizing 
the effects of the interaction in the CP equation is given by dN. This implies that 
one can get the same results with a proper rcscaling of the variables N and a. As 
it can be seen by a simple inspection of the equations, the scaling property is lost 



Table 1 



Chemical potential /x, and ground state energy per particle e, of Rb atoms in an 
isotropic trap {uj/2'k = 77.78Hz) in different approaches. Energies are in units of Tuo. 
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TF GP MGP CHhnc 


TF GP MGP GHhxc 


10^ 
10^ 
10^ 


16.75 16.85 16.99 16.94 
42.07 42.12 42.69 42.53 
105.68 105.70 107.97 107.20 


11.96 12.10 12.19 12.20 

30.05 30.12 30.48 30.48 
75.49 75.52 76.94 76.85 



in the MGP approach. Moreover, the relation between the different contributions to 
the chemical potential changes to 



^l = ekin + eno + '^e^l + ^e\2l, 



,(2) 
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where 
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In this case the relation implied by the virial theorem is 

2eki„ - 2eHO + ^€1 + ^€1 = 0- 



*i(f). 



(17) 
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A simple approach, valid for large N and loosely called the Thomas-Fermi (TF) 
approximation, is obtained by neglecting the kinetic energy term in the GP equation. 
In the TF approximation it is possible to derive simple analytical expressions [14] , use- 
ful to make quick estimates of several quantities. For instance, //tf = l/2{15dN)'^/^, 
while the energy is related to the chemical potential by btf = 5/iTF/7. The local 
value of a; = Na^pi{0) at the center of the density distribution of the trapped bosons 

is given by a;TF(0) = /{Stt). 

As can be seen in Tabic I, for large number of particles the TF and GP re- 
sults are practically identical and there is also a very good agreement between the 
MGP and GHhnc results. In the case = 10^, the contributions to the energy 
for the GP(MGP) equations are: eki„ = 0.0294(0.0292), cho = 45.306(46.57) , 
e^l = 30.184(28.96) and Cint = 0(1.379). The virial theorem is well fulfilled in 
both cases . 

By changing the number of trapped atoms in the range of the experimental 
availability, the average values of x are such that the corrections to the GP equation 
are kept small and of the order of 2% in the case of = 10^ atoms. However, the 
recent experiments, where the scattering length can be largely manipulated, open 



the door to explore higher values of x. In fact, in order to vary x, it is much more 
efficient to change the scattering length than the number of atoms. Experimental 
results are available for ^^Rb, whose scattering length is modulated by its Feshbach 
resonance. The number of trapped atoms is N ^ 10^ and the trap is anisotropic. In 
order to estimate the corrections induced by the MGP equation, we have considered 
an isotropic trap characterized by the frequency u)/2tt = 10 Hz, corresponding to an 
average of the cylindrical trap frequencies used in the experiment. The u value is 
smaller than that employed for ^^Rb and therefore oho is larger. We take a — 0.1228, 
which is in the range considered by the experiments. It corresponds to a = SOOOao, 
where oq is the Bohr radius of the hydrogen atom. In this case, xtf(O) ~ 0.03, which 
is just beyond the range of the points plotted in Fig.l. The energies per atom turn 
out to be: cqp = 18.25 and cmgp — 21.85. For the chemical potential we have: 
piQ-p = 25.48 and /xmgp = 31.09. As a consequence of the use of the MGP equation, 
the corrections are of the order of 20 %. The actual cylindrical trap is currently 
under study [15]. 



4. COLLECTIVE EXCITATIONS 

Information on the excitation spectrum of a system are contained in the dynamic 
structure function, which, for a given excitation operator F, is given by 

Sf{E) = |(n|F|0)p5(£; - {En - Eo)). (19) 

n 

In the case of inelastic neutron scattering ( for instance against liquid ^He) the 
operator F corresponds to the density fluctuation operator, F — e^'^^J , and the 
inclusive inelastic cross section is proportional to the dynamic structure function. 

A useful tool to anlyze Sf{E) is provided by the sum rules approach, extensively 
used in quantum many-body systems both in the context of nuclear [16] and con- 
densed matter physics, in particular to analyze the excitation spectrum of quantum 
liquids [5]. The sum rules establish rigorous links among the energy momenta of 
Sf{E) and ground state properties. The energy weighted integrals are defined as: 

ruk^ E^SF{E)dE, (20) 
Jo 

and ruk can be calculated, without an explicit knowledge of Sf{E), as a ground 
state expectation value of certain combinations of commutators of the excitation 
operator F and the Hamiltonian. So, quantities like rrik+i/mk or (771^+2/""^/=)^^^ 
supply rigorous upper bounds to the energy of the lowest excited state that can be 
connected to the ground state through the operator F. The upper bounds are very 
reliable when the excited state is highly collective, i.e., when the strength distribution 
is almost exhausted by a single mode. We will concentrate here in the case of the 
compressional modes, or monopole excitations, whose associated excitation operator 

• 7-1 v^A^ 2 

isF = ^^ rf. 

In the present case, we will study the upper bound provided by (ma/mi)^/^. 
Using the definition of the dynamic structure function and the completness of the 



eigenstates \n), the moments mi and 777,3 can be expressed as 



mi = ^(0|[Ft,[iy,F]]|0) 



(21) 



and 



m, = l{0\[[F\H],[H,[H,F]m. 



(22) 



By explicitly calculating the commutators, mi is expressed in terms of the mean 
square radius, 
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mi = 2— 7V(0|r^|0) . 
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(23) 



Obviously, we get the same expression of nii for both GP and MGP equations. 
In the case of the monopole excitation, is more efficient to calculate ms as 
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2 V m j dX^ 



(24) 



where A is the parameter of the scaling transformation p(r) X^p{Xr). In fact, for 
the GP equation, the energy for the scaled density is 



E{p, A) = iV [x'ekM + ^ + >^'4nt 



(25) 



The condition dE{p, \)/dX = at A = 1 satisfies the virial theorem, ms is given by 
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Finally 
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By a similar procedure for the MGP equation we get 



e^- 27 e^'^ 



+ 



eHO 8 cho 



1/2 



(28) 



In the ^^Rb case the estimates of the monopole excitation energies (in HO units) 
are 2.23 and 2.38 for the GP and MGP, respectively. So, the MGP correction to 
the excitation energy is about 7%. This correction is smaller than that to the energy 
itself. The whole spectrum is shifted but the separation between the excitation energy 
levels is less affected. 



4. CONCLUSIONS 



In conclusion, we find that the MGP equation induces corrections of 20% in the 
ground state properties of the condensate, when the conditions of the recent exper- 
iments for ^^Rb are considered [x ~ 10~^). MGP is still a mean field theory, since 
it tries to incorporate correlation effects into the average single particle potential, 
and it cannot predict the depletion of the condensate. However, we believe that the 
estimates of the energy, chemical potential, and density profile are surely indicative 
and can still be accurate. Moreover, it is legitimate at these densities to question 
the use of a simplified interaction, given in terms of hard spheres. In any case, it is 
clear that fully microscopic calculations, which might address the many-body wave 
function and take explicitly into account the depletion of the condensate, are urgently 
required [17]. 
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